The disclosure of the invention includes a microfiche attachment containing source code for computer programs described herein. Microfiche Attachment A, which consists of three fiche (Nos. 1-3) with a total of 153 frames, is incorporated herein by reference in its entirety. The Computer Appendix contains material that is subject to copyright protection. The copyright owner has no objection to the facsimile reproduction by anyone of the patent document or patent disclosure, as it appears in the Patent and Trademark Office patent file or records, but otherwise reserves all copyright rights.
The present invention relates generally to signal encoding and reconstruction. More particularly, the present invention relates to devices and methods for reconstructing a signal from measurement data from a detector.
Exact measurement of the properties of nature is a common goal within the experimental sciences. Similarly, medical diagnostics and communications technology, among other scientific endeavors, seek the ability to obtain exact measurement of properties within their respective fields, e.g., MRI or free space optical or other electromagnetic transmission. Optimal extraction of the underlying measurement data requires the removal of measurement defects such as noise and limited instrumental resolution. However, in spite of the availability of highly sophisticated instruments, noise and instrumental signatures are still present in the data, making the measurement only approximate.
An area of experimental science in which instrumental signatures are particularly a problem is astronomy, where the signal sources to be measured are very faint. Even when the instruments are made essentially noise-free, instrumental signatures related to finite spatial, spectral, or temporal resolution remain. At this point, image reconstruction is required to remove the instrumental signatures.
One of the most powerful approaches to image restoration, e.g., removal of point-spread-function blurring, is Bayesian image reconstruction, which includes Goodness-of-Fit (Maximum Likelihood) and Maximum Entropy data fitting. This family of techniques employs a statistical relationship between various quantities involved in the imaging process. Specifically, the data, D, consisting of the original noisy, blurred image model is linked to the desired noise-free, unblurred image model, I, through a model, M. The model M includes all aspects of the relationship between the data and the image, e.g., that the data is normally collected on a rectangular grid and that the data is related to the image model through the relationship
D(i)=(*H)(i)+N(i)xe2x80x83xe2x80x83(1)
where D(i) is the data in cell i (typically a pixel), I is the image model, H is the point-spread-function (PSF) due to instrumental and possible atmospheric blurring, * is the spatial convolution operator, i.e.,                                           (                          f              *              g                        )                    ⁢                      (            x            )                          =                              ∫                          -              ∞                        ∞                    ⁢                      xe2x80x83                    ⁢                                    ⅆ                              x                xe2x80x2                                      ⁢                          f              ⁡                              (                                  x                  xe2x80x2                                )                                      ⁢                          g              ⁡                              (                                  x                  -                                      x                    xe2x80x2                                                  )                                                                        (        2        )            
and N represents the noise in the data, assuming the PSF is a function only of displacement between pixels. In general, the PSF can vary across the field.
Image reconstruction differs from standard solutions of integral equations due to the noise term, N, that nature of which is only known statistically. Methods for solving such an equation fall under the categories of (1) direct methods, which apply explicit operators to data to provide estimates of the image, and (2) indirect methods, which model the noiseless image, transform it only forward to provide a noise-free data model, then fit the parameters of the image to minimize the residuals between the real data and the noise-free data model. The direct methods have the advantage of speed, but they tend to amplify noise, particularly at high frequencies. The indirect methods supposedly exclude the noise, however, the required modeling is a disadvantage. If a good parametric form for the image is known a priori, the result can be very good.
To statistically model the imaging process, Bayesian image reconstruction methods analyze the properties of the joint probability distribution of the triplet, D, I and M, i.e., p(D,I,M). Applying Bayes"" Theorem [p(A,B)=p(A B)p(B)=p(B A)p(A), where p(X Y) is the probability of X given that Y is known] provides:                                                                         p                ⁡                                  (                                      D                    ,                    I                    ,                    M                                    )                                            =                                                p                  ⁡                                      (                                                                  D                        ⁢                                                  xe2x80x83                                                ⁢                        I                                            ,                      M                                        )                                                  ⁢                                  p                  ⁡                                      (                                          I                      ⁢                                              xe2x80x83                                            ⁢                      M                                        )                                                  ⁢                                  p                  ⁡                                      (                    M                    )                                                                                                                          =                                                p                  ⁡                                      (                                                                  I                        ⁢                                                  xe2x80x83                                                ⁢                        D                                            ,                      M                                        )                                                  ⁢                                  p                  ⁡                                      (                                          D                      ⁢                                              xe2x80x83                                            ⁢                      M                                        )                                                  ⁢                                  p                  ⁡                                      (                    M                    )                                                                                                          (        3        )            
By setting the first factorization of p(D,I,M) in Equation 3 equal to the second factorization provides the usual starting point for Bayesian reconstruction:                               p          ⁡                      (                                          I                ⁢                                  xe2x80x83                                ⁢                D                            ,              M                        )                          =                                            p              ⁡                              (                                                      D                    ⁢                                          xe2x80x83                                        ⁢                    I                                    ,                  M                                )                                      ⁢                          p              ⁡                              (                                  I                  ⁢                                      xe2x80x83                                    ⁢                  M                                )                                                          p            ⁡                          (                              D                ⁢                                  xe2x80x83                                ⁢                M                            )                                                          (        4        )            
A common goal of Bayesian image reconstruction is to find the M.A.P. (Maximum A Posteriori) image, I, which maximizes p(I D,M), i.e., the most probable image given the data and model. (Note that other image estimates, e.g., the average image,  less than I greater than =∫D,M dMdD I p(I D,M), may be used here and in the methods described in the detailed description.) (MAP image reconstruction is also sometimes known as parametric least-squares fitting.)
It is common in Bayesian image reconstruction to assume that the model is fixed. In this case, p(D M) is constant, so that
p(I D, M)∞p(D I,M)p(I M).xe2x80x83xe2x80x83(5)
The first term, p(D I,M), is a goodness-of-fit quantity, measuring the likelihood of the data given a particular image and model. The second term, p(I M), is normally referred to as the xe2x80x9cimage priorxe2x80x9d, and expresses the a priori probability of a particular realization of the image given the model. In Goodness-of-Fit (GOF) image reconstruction, p(I M) is effectively set to unity, i.e., there is no prior bias concerning the image. Only the Goodness-of-Fit (p(D I,M)) is maximized during image reconstruction. Typically,
p(I|D,M)=exp (xe2x88x92X2R/2),xe2x80x83xe2x80x83(6)
where X2R is the chi-square of the residuals, R (≈Dxe2x88x92I*H). While this approach ensures that the frequency distribution of the residuals has a width which is characteristic of the noise distribution, it normally results in images with spurious spatial features where the data has a low signal to noise ratio (SNR). Also, the large amplitude residuals often show a strong spatial correlation with bright features in the data.
When no parametric model of the image is known, the number of image model parameters can quickly become comparable to, or exceed, the number of data points. In this case, a MAP solution becomes inappropriate.
For example, if the number of points in the image model equals the number of data points, the nonsingular nature of the transform in Equation 1 assures that there is a solution for which the data, including the noise, are exactly modeled with zero residuals. This is the same poor solution, with all its noise amplification, obtained by the naive Fourier deconvolution. Thus, an unrestricted indirect method is no better at controlling noise than a direct method, and therefore, the image model must be restricted in some way. The indirect methods restrict the image model and differ only in the specifics of image restriction.
A simple restriction is to constrain the model image to be positive. Since even a delta-function image is broadened by the PSF, it follows that the exact inverse of any noisy data with fluctuations on scales smaller than the PSF must be both positive and negative. By preventing the image from becoming negative, the noise-free data model cannot fluctuate on scales smaller than the PSF, which is equivalent to smoothing the data on the scale of the PSF. However, this Non-Negative Least-Squares (NNLS) fit method is not able to eliminate noise fitting on larger scales.
Maximum entropy (ME) image reconstruction solves many of the problems of the simpler GOF methods (e.g., NNLS). In ME imaging, one calculates a value for the image prior based upon xe2x80x9cphase space volumexe2x80x9d or counting arguments. Heuristically, p(I M) is written p(I M)=exp(S), where S is the entropy of the image in a given model. All ME methods capitalize on the virtues of incorporating prior knowledge of the likelihood of the image. The information entropy is maximized for a flat image, which eliminates structure not required by the data and suppresses noise fitting. The benefits of this are numerous, including eliminating the over-resolution problems of GOF methods and increasing the numerical stability of the calculations. The disadvantages are that the image prior used in ME takes the form of a global constraint, resulting in the enforcement of an average smoothness (average information content) on the entire image without recognizing that the density of information content, and, thus, the appropriate degree of smoothing, varies from location to location in the image. This results in oversmoothing of the image in some locations and undersmoothing it in others.
Many Bayesian image reconstruction methods assume that the model is fixed. Others, however, propose varying the model, such as the multi-channel image modeling of Weir (Applications of Maximum Entropy Techniques to HST Data, Proceedings of the ESO/ST-ECF Data Analysis Workshop, April 1991). In this method, the image is assumed to be a sum of pseudoimages convolved with a blurring function of various spatial scales. This method, while superior to many of its predecessors, may exhibit low- level spurious sources as ripples in the reconstructed image, and still displays some spatial correlation within the residuals.
The need remains for a method of image reconstruction which is capable of effectively reducing noise fitting without relying on a priori probabilities that result in excessive averaging of the image, but rather which adapts itself to the distribution of information content in the image. The following disclosure provides such a method.
It is an advantage of the present invention to provide a method for identifying a generalized image cell as an optimal basis for image reconstruction.
It is another advantage of the present invention to provide a method for minimizing the number of parameters, i.e., minimizing image model complexity, needed to reconstruct an image.
It is still yet another advantage of the present invention to provide a method for accelerating the computational speed for image reconstruction.
The method identifies a Pixon element, which is a fundamental and indivisible unit of information, and a Pixon basis, which is the set of possible functions from which the Pixon elements are selected. The actual Pixon elements selected from this basis during the reconstruction process represents the smallest number of such units required to fit the data and represent the minimum number of parameters necessary to specify the image. The Pixon kernels can have arbitrary properties (e.g. shape, size, and/or position) as needed to best fit the data. The only criterion is that the kernel functions be selected to efficiently model the image with the fewest number of components. The goodness-of-fit (GOF) is determined according to the chi-squared statistic, the Maximum Residual Likelihood (MRL) statistic, or a similar chi-squared statistic as appropriate to the data set being modeled.
The goal of the Pixon method of image reconstruction is to find the simplest model meeting a given quality of fit using the chosen GOF criterion. Unlike standard Bayesian methods, the Pixon method need not explicitly define an image prior (although it could). Rather it recognizes that finding the minimum complexity model for a given quality of fit satisfies the general goal of Bayesian methods, since simple models give rise to optimized image priors. If a prior were to be defined, it would increase monotonically with image model complexity. In this case, the specific functional form for the prior would have to be selected for the specific imaging case under consideration, e.g. based on the physics of the situation or built up from statistics gathered from a training set of data. Normally, however, it is sufficient to have a criterion with which to decide which is the simplest model among a family of models and to use this to choose the simplest model from amongst all images fitting the data to the specified quality.
The benefits of building a minimum complexity model are many. Since a minimum complexity model more critically fits the image to the data, the parameters of the image are more accurately determined since a larger fraction of the data is used to determine each one. For the same reason, a minimum complexity model does not show signal-correlated residuals, and hence provides unbiased source strength measurements to a precision limited only by the theoretical limits set by the noise statistics of the data. In addition, since the image is constructed from a minimum complexity model, spurious (i.e., artificial or numerically created) sources are eliminated. This is because a minimum complexity model only has sufficient parameters to describe the structures that are required by the data and has none left over with which to create false sources. Finally, because the Pixon method builds a critical model and eliminates spurious sources, it can achieve greater spatial resolution than competing methods and detect fainter sources that would otherwise be hidden by these spurious structures.
In application to astronomy, the input data is commonly obtained with a CCD detector array with a particular pixel size and shape. For example, in the imaging of star fields, the image would be best represented as a sum of point sources with arbitrarily precise positions and brightness. Since large regions of the data field will have few, if any, photon counts, portions of the pixel grid are unused, and the degrees of freedom, i.e., pixels, representing these portions of the image over-specify the data. In other portions of the image, the density of pixels may be too sparse to adequately represent the image. In this example the sizes, shapes of the Pixon kernels would be adjusted to decompose the image into the fewest number of Pixon elements. In Bayesian terms, since the number of Pixon kernels used would be much fewer than the total number of pixels in the grid, any sensible image prior term in Equation 5 (in the Background of the Invention) would be greatly improved.